OPTIMALLY SWIMMING STOKESIAN ROBOTS 
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LEFEBVRE-LEPOT§ , AND BENOiT MERLET^ 

Abstract. We study self-propelled stokesian robots composed of assemblies of balls, in dimen- 
sions 2 and 3, and prove that they are able to control their position and orientation. This is a result 
of controllability, and its proof relies on applying Chow's theorem in an analytic framework, similar 
to what has been done in [3] for an axisymmetric system swimming along the axis of symmetry. We 
generalize the analyticity result given in [3] to the situation where the swimmers can move either in 
a plane or in three-dimensional space, hence experiencing also rotations. We then focus our atten- 
^~v tion on energetically optimal strokes, which we are able to compute numerically. Some examples of 
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computed optimal strokes are discussed in detail 
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Q 1. Introduction. Self-propulsion at low Reynolds number is a problem of con- 

siderable biological and biomedical relevance which has also great appeal from the 
point of view of fundamental science. Starting from the pioneering work by Taylor 
'""' [32] and Lighthill [23], it has received a lot of attention in recent years (see the en- 

I— I cyclopedia article [4] for an elementary introduction and the review paper |22j for a 

Q_J comprehensive list of references) . 

Both relevance for applications and theoretical interest stem from the fact that 
one considers swimmers of small size. The Reynolds number Re = LV/v gives an es- 
timate for the relative importance of intertial to viscous forces for an object of size L 
2 moving at speed V through a newtonian fluid with kinematic viscosity v. Since in ap- 

^ plications V rarely exceeds a few body lengths per second, if one considers swimming 

in a given medium, say, water, then Re is entirely controlled by L. At small L, inertial 
fvj forces are negligible and, in order to move, micro-swimmers can only exploit the vis- 

^ cous resistance of the surrounding fluid. The subtle consequences of this fact (which 

^^ are rather paradoxical when compared to the intuition we can gain from our own 

vN swimming experience) are discussed in |30j . For example, the motion of microswim- 

mers is geometric: the trajectory of a low Re swimmer is entirely determined by the 
sequence of shapes that the swimmer assumes. Doubling the rate of shape changes 
1^- simply doubles the speed at which the same trajectory is traversed. As observed in 

|31j . this suggests that there must be a natural, attractive mathematical framework 
for this problem (which the authors, indeed, unveil). On the other hand, bacteria 
?• and unicellular organisms are of micron size, while artificial robots to be used non- 

. ^H invasively inside human bodies for medical purposes must be small. Discovering the 

j^ secrets of biological micro-swimmers and controlling engineered micro-robots requires 

J-^ a quantitative understanding of low Re self-propulsion. 

The basic problem of swimming is easy to state: given a (periodic) time history 
of shapes of a swimmer (a sequence of strokes), determine the corresponding time 



*CMAP UMR 7641, Ecolc Polytochnique CNRS, Route de Saclay, 91128 Palaiseau Cedex - France 
f rancois . alougesapolytechnique . edu) 

^SISSA, International School of Advanced Studies, Via Bonomea 265, 34136 Trieste - Italy 
desimoneSsissa. it) 

■•"SISSA, International School of Advanced Studies, Via Bonomea 265, 34136 Trieste - Italy 
luca.heltaiOsissa. it) 

§CMAP UMR 7641, Ecole Polytechnique CNRS, Route de Saclay, 91128 Palaiseau Cedex - France 
aline . lef ebvreOpolytechnique . edu) 

IICMAP UMR 7641, Ecolc Polytochnique CNRS, Route de Saclay, 91128 Palaiseau Cedex - France 
benoit .merletOpolytechnique .edu) 



history of positions and orientations in space. A natural, related question is the 
following: starting from a given position and orientation, can the swimmer achieve 
any prescribed position and orientation by performing a suitable sequence of strokes? 
This is a question of controllability. The peculiarity of low Re swimming is that, since 
inertia is negligible, reciprocal shape changes lead to no net motion, so the question 
of controllability may become non trivial for swimmers that have only a few degrees 
of freedom at their disposal to vary their shape. The well established scallop theorem 
[50] is precisely a result of non-controllability. 

Once controllability is known, i.e., it is shown that it is possible to go from A to 
B, one can ask the question of how to go from A to B at minimal energetic cost. This 
is a question of optimal control. 

In spite of the clear connections between low Re self-propulsion and control theory, 
this viewpoint has started to emerge only recently, and mostly in the mathematical 
literature. Examples are [51], and the more recent contributions [TU], [H], [20] and 
[26] ■ The papers [1], [3], [1] study in detail both controllability and optimal control 
for axisymmetric swimmers whose varying shapes are described by few (in particular, 
two) scalar parameters. 

In this paper we analyze the problem of low Reynolds number swimming from the 
point of view of geometric control theory, and focus on a special class of model systems: 
those obtained as assemblies of a small number of balls. Model systems of this type 
have played an important role in clarifying the subtleties of low Re self-propulsion 
[27] ■ [6], [13]. These systems offer an interesting balance between complexity of the 
analysis and richness of observable behavior. 

While constructing such artificial swimmers may indeed be possible, their prac- 
tical uses are at present unclear. However, the analysis of their swimming patterns 
may prove a very useful tool both in the design of robotic microswimmers |14| , and in 
understanding the motion of biological swimmers. As an example, the social motility 
patterns exhibited by Myxobacteria (see, e.g., [18]) are strikingly similar to the ones 
described in this paper. As part of their life cycle, individual cells become linked 
together and move collectively on surfaces. The links between adjacent cells are guar- 
anteed by pili which individual cells can project and retract. The fact that collective 
motion rests upon control of the relative positions of individuals seems fully estab- 
lished. Our study may shed light on the many details that are at present unclear. 

Proving controllability and providing optimal control strategies for realistic swim- 
mers, either biological or artificial, is rather difficult. The use of assembly of balls as 
a model swimmer relies on the idea that one can use collections of spheres to model 
the hydrodynamic characteristics of real swimmers, replacing their actual swimming 
apparatus (which can be indeed very complex to describe, simulate and approximate) 
with a finite collection of spheres having comparable hydrodynamic resistance. Very 
recent progress towards the analysis of low Re more general swimmers propelling 
through shape changes can be found in [T^] , [IS] . 

Our approach is similar in spirit to the one in [5], [3]j [l]j but we extend it to 
non-axisymmetric systems such as three spheres moving in a plane, and systems of 
four spheres moving in three dimensional space. The motion of these systems is 
described by both positional and orientational variables, leading to a much richer 
geometric structure of the state space. The study of such systems requires substantial 
extensions of the mathematical and numerical methods introduced in [3] . 

For all the model swimmers described above, controllability is proved by using 
two main ingredients. The first is Chow's theorem, leading to local controllability in 
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a neighborhood of a point X in state space. We verify the full rank (5.2) hypothesis 
of this theorem by showing that the vector fields of the coefficients of the governing 
ODEs (a linear control system without drift) and their first order Lie brackets span 
the whole tangent space to the state space at X. The second is the analyticity of the 
coefficients, and the fact that our shape space is always connected. This allows us 
to pass from local to global controllability generalizing a result that has been proved 
in [3] for the special axisymmetric case. In the same paper [35, Chow's theorem was 
also used to prove local controllability. In this simpler axisymmetric case, however, 
position is described only by one scalar parameter (the position p of one distinguished 
point along the axis of symmetry). The non-degeneracy condition reduces then to 
the non-vanishing of a single scalar quantity, namely, the curl of the vector field V, 
governing the rate of change of position as a consequence of shape changes at rates 
^i according to p = Vi(^)^i + V2(^)^2- In the richer context of this paper, which 
involves also rotational degrees of freedom of the swimmers, proving controllability 
requires an explicit computation of all the first order Lie Brackets. In fact, all the 
systems we analyze here satisfy the condition 

M + M(Af-l)/2 = d 

where M is the number of controls (rate of change of shape variables) and d is the 
dimension of the state space (position, orientation, and shape). This is the necessary 
condition that first order Lie brackets alone suffice to show that the Lie algebra of the 
coefficients has full rank, so that controllability follows. 

For controllable systems, it makes sense to ask how to achieve the desired target 
(position and orientation) at minimal energy cost. We present a method to address 
this optimal control question numerically, and we then examine in detail several op- 
timal strokes for a concrete model swimmer (three balls swimming in a plane with a 
prescribed lateral displacement). Depending on whether the final orientation is also 
prescribed, and whether the initial shape is prescribed as well or rather one treats 
it as a parameter to be optimized, we obtain dramatically different answers. Their 
variety illustrates the surprisingly richness of behavior of low Re swimmers. 

The rest of the paper is organized as follows. In Section[2] we describe the various 
model swimmers to which our analytical and numerical tools are later applied. The 
first is the Najafi-Golestanian's swimmer P^, already treated in [3], while the two 
others are non-trivial generalizations. Section [3] presents some results on Stokes flows 
and in Section |4] we show that swimming is indeed an affine control problem without 
drift. In SectionlSJ we prove the effective swimming capability of our model swimmers. 
In Sections [6] and [7] we state the optimal control problem and a numerical strategy 
for its solution. Examples of optimal strokes for three balls swimming in a plane are 
discussed in detail in Section [H 

2. The svifimmers. We will focus our attention on some special swimmers. 
Namely, we assume that the swimmer is composed of A^ non-intersecting balls (-Bi)i<i<Ar 
C M? centered at (xi)i<j<jv, and we restrict ourselves to configurations which can be 
described by two sets of variables: 

• the shape variables^ denoted by ^ S 5, where S is an open connected subset of 
R*^, from which relative distances {xij)i<:ij<:pf between the balls (-Bi)i<i<Ar 
are obtained. In the examples treated in this paper, the balls are assumed 
to move only along fixed directions which make fixed angles one to another. 
This reflects a situation where the balls are linked together by thin jacks 
that are able to elongate. The viscous resistance associated with these jacks 
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Fig. 2.1. Three Sphere Swimmer (3S) 



is, however, neglected, and the fluid is thus assumed to fill the whole set 

• the position variables, denoted hy p £ V, which describe the global position 
and orientation in space of the swimmer. In our examples, the set V is typi- 
cally a manifold of dimension less than or equal to six. The six-dimensional 
case is given hy p G V = R^ x S0(3) and p consists of a translation and a 
rotation in the three-dimensional space. 
We also assume that the orientation of the balls (i3i)i<i<Ar and the distances 
{y:.ij)i<ij<N depend analytically on (^,p), therefore the state of the system is analyt- 
ically and uniquely determined by the variables (^,p), so that there exist TV (analytic) 
functions 

X, : S X V X dB — ^ R3 

which give the position of the current point of the i— th sphere of the swimmer in 
the state (^,p) in 5 x 7^. The non-slip boundary condition on dBi imposes that the 
velocity of the fluid is given by 

v,(C,p,r) = ^X,(C,p,r) - ii ■ V5)X,;(e,p,r) + (p ■ Vp)X,(^,p, r). (2.1) 

2.1. The three sphere swimmer of Najafi and Golestanian (3S). This 
swimmer, initially proposed in |27j has been studied thoroughly in j3J and j5]- It is 
composed of 3 spheres of radius a > aligned along the x— axis, as depicted in Fig. 

We call S, = (^1,^2) the length of the arms, and p the position of the central 
sphere which leads us, in order to avoid overlap of the spheres, to consider the shape 
set S = (2a, -l-oo)^ and the position set P = R. Indeed, due to axial symmetry, this 
swimmer may only move along the x— axis. Using ei = (1, 0, 0)"^, we can write 

xi(C,p) = (p-6)ei,X2(C,p) =pei ,X3(^,p) = {p + £,2)ei, 

and 

X,(C,p,r)=x,(C,p) + r, Vte {1,2,3}, VredB. 

2.2. The three sphere swimmer moving in a plane (3SP). A variant of 
the Najafi-Golestanian's swimmer, called "Purcell's rotator", was presented in [13], 
where the axis along which the three spheres are constrained to move, is bent to 
form a circle with fixed radius. The resulting swimmer is one where the three spheres 
always keep a fixed distance from the center of the swimmer, and may only vary their 
location on the circle, remaining on the same plane. 
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Fig. 2.2. The three-sphere planar swimmer (3SP) 



As observed in [13_, this variation allows one to control rotations around the axis 
perpendicular to the circle and passing through the center of the swimmer (hence the 
name ^^rotator™), and introduces as well a small drift in the system, which displaces 
horizontally on the plane containing the circle. 

That swimmer does not change the number of controls (two) with respect to 
the Najafi Golestanian's swimmer, however a periodic shape change induces both a 
rotation and a translation of the swimmer (that is, the dimension d of the state space 
is equal to 5 ^ Af + M{M — l)/2 = 3). Seen from the perspective of control theory, 
this implies that not only first order Lie brackets are used to change the positional 
variables but also higher order ones, with smaller and smaller effectiveness as the 
order of the used Lie brackets increases (see Section [s]) . 

A complementary version of Purcell's rotator, with one additional control variable, 
has been proposed and studied in [53] ■ It is composed of three balls Bi, B2, B3 of equal 

radii a > 0. The three balls c an m ove along three horizontal axes that mutually meet 

2tt 



at c with angle — - (see Fig. 



2.21). The balls do not rotate around their axes so that 

the shape of the swimmer is characterized by the three lengths ^1,^2,^3 of its arms, 
measured from the origin to the center of each ball. However, the swimmer may freely 
rotate around c in the horizontal plane. 

Consider a reference equilateral triangle {81,82, S3) with center O G R'^ in the 
horizontal plane {0,x,y) such that dist(0,S'i) = 1 and define tj = 08i. Position 
and orientation in the horizontal plane are described by the coordinates of the center 
c e R'^ (but c stays confined to the horizontal plane) and the horizontal angle 9 that 
one arm, say arm number 1, makes with a fixed direction, say {0,x), in such a way 
that d — i. Therefore, we place the center of the ball Bi at x^ = c + S^iTZeti with 
^i > for i = 1,2, 3, where TZg stands for the horizontal rotation of angle 9 given for 
instance by the matrix: 

cos(6i) - sin(6') 

7^fl = I sin(6') cos{9) 

1 



The swimmer is then fully described by the parameters X 
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i^,c,9) eSxV, 



where S :— (-%, +oo)'^, the lower bound being chosen in order to avoid overlaps of 
the balls, V = R^ x R, and the functions X,; are now defined as 

X,(e, c, a, r) = c + 7^e(6t, + r) Vi e {1, 2, 3} . 

Notice that the functions X^ are still analytic in (^, c, 9), and we use them to compute 
the instantaneous velocity on the sphere Bi 

9X,; 



dt 



■ (C, c,e,r) = c + Oes x (^.t, + r) + UeU^r , 



where 63 is the vertical unit vector. Eventually, due to the symmetries of the system, 
the swimmer stays in the horizontal plane. 

2.3. The four sphere swimmer moving in space (4S). We now turn to 
the more difficult situation of a swimmer able to move in the whole three dimensional 
space and rotate in any direction. In this case, we fix TV = 4 and we consider a regular 
reference tetrahedron (S"!, 5*2, S'3, 5*4) with center O G R^ such that dist(0,S'i) = I 
and as before, we call t^ = OSj for i = 1,2, 3, 4. 

The position and orientation in the three dimensional space of the tetrahedron 
are described by the coordinates of the center c g R'^ and a rotation 7?. G SO (3), in 
such a way that d ~ 6. 

We place the center of the ball i?i at x^ = c + ^{R-ti with ^^ > for z = I, 2, 3, 4 
as depicted in Fig. |2.3| and forbid possible rotation of the spheres around the axes. A 
global rotation [TZ 7^ Id) of the swimmer is however allowed. 

The four ball cluster is now completely described by the list of parameters X = 

(C,c,7^) G S xP, where S := (Vf , +00)'' and P = R^ x S0(3). Again, the lower 
bound for ^^ is chosen in order to avoid overlaps of the balls. 
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Fig. 2.3. The four sphere swimmer (4S). 

Furthermore, the function X^ are now defined as 

X,(^, c, 7^, r) = c + n{i,t, + r) Vz e {1, 2, 3, 4} , 

which are still analytic in (^, c, TV), from which we compute the instantaneous velocity 
on the sphere Bi 



dt 



- (e, c, 7^, r) = c + u; X (^.t, + r) + UtS 



where a; is the angular velocity. Recall that u) — Ax{lZTZ) is the axial vector associated 
with the skew matrix TZTZ. 



2.4. The result. In the theoretical part of this paper, we will establish the 
following theorem. 



Theorem 2.1. Consider any of the swimmers described in Sections 2.1 2.2 
or \2.3l and assume it is self-propelled in a three dimensional infinite viscous flow 
modeled by Stokes equations. Then for any initial configuration (^',p*) € S x V any 
final configuration {£,■' ,p-^) ^ S x V and any final time T > 0, there exists a stroke 
C G C°{[0,T],S), piecewise Ci([0,T],5) satisfying ^(0) = C and ^{T) = ^^ such that 
if the self-propelled swimmer starts in position p^ with the shape ^* at time t — 0, it 
ends at position pf and shape S,^ at time t — T by changing its shape along £,{t). 

3. Modelization of the fluid. In this section we give the expression of the 
total force and torque induced by a shape change of the swimmer. Since the swimmer 
is composed of unions of balls, this turns out to study the Dirichlet-to-Neumann map 
outside the swimmer. We recall that the swimmer is composed of N identical and 
non-intersecting balls Bi , • • • , Bn of radius a > not necessarily aligned, of center 
Xi e R'^, but linked together by deformable jacks which form a kind of skeleton. Since 
the lateral size of the arms of this skeleton are negligible, we can consider that the 
fluid fills the unbounded domain ft = R'^ \ ufLj^Bi. Since the balls do not intersect, 
the vector x = (xi, • • • , xjv) € fi is restricted to belong to 

Sn ■■= |xe(R3)^ : min|x, -Xjl > 2ai . 

We work at low Reynolds number, so that the fluid obeys Stokes equations 
-77AU + Vp == mft, 

^■" = ? '"^o (3-1) 

—an = t on oil, 

u ^^ at 00 . 

where (u,p) are respectively the velocity and the pressure of the fluid, rj its viscosity, 
(T := ?7(Vu + Vu*) — pid is the Cauchy stress tensor and n is the outer unit normal 
to dil (hence, n points from the fluid to the interior of the balls). Existence and 



uniqueness of a solution to (3.1 1 is classical in the Hilbert space 

u 



V := \ue V'{n, r3) I Vu e L^{n), e L^n) \ , 



endowed with the norm 



iVup 



Assuming that the force field belongs to Ti ^/^(9rirl the solution of (3.1) can be 



expressed in terms of the associated Green's function, namely the stokeslet 



1 /I r ® r^ 



G(r) := ^(0+^^ (3-2) 



^Here, and in all the paper the symbol H" denotes the Sobolev space of order s. 
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as 



u(r) = f G(r-cr)f(cr)dcr 
Jan 



where f is a distribution of forces on dfl. The Neumann-to-Dirichlet map 

f ^ U|ao= { I G{v-a)i{a)da] 

\Jdn J |an 



(3.3) 



(3.4) 



is a one to one mapping onto while, by the open mapping Theorem, its inverse (the 
Dirichlet-to Neumann map) 7~~^ is continuous. 

Calhng B = B{0,a) the ball of radius a centered in the origin O, the special 
structure of our domain allows us to identify dfl with [dB)^ in such a way that T 
can be seen as a map 

/x -H^ -^ Hn (3.5) 

(fi,--- ,f7v) I — > (ui,--- ,UAr) 

in which H^ stands for {H''{dB)) and the dependence on x S Sn has been empha- 
sized, defined by 

N 

Uj(r) = y] / G(xi -Xj +r -o-)fj(cr)d(T 

^=1 (3.6) 

=: ^(f„G(x,:-x,+r-.))gB' ^ € 95. 

3.1. Analyticity of the Dirichlet-to-Neumann map. This section is de- 
voted to the following lemma which will be essential in the proof of the controllability 
theorem. This kind of result was already proved in [3] but with a much more compli- 
cated method. We therefore give another much simpler proof which applies to more 
general situations, e.g., the swimmers described in subsections |2.2| and |2.3[ We de- 
note by C{E, F) the Banach space of linear maps from E to F endowed with its usual 
norm. 

Lemma 3.1. The mapping x i — > Tx is analytic from Sn into C{Hj^ iT~^n )■ 

The notion of analyticity in a Banach space is classical: it means that at all points 
Xq S Sn^ Tx is equal to its Taylor series which converges in L{'H~^ i^jv ) ^^^ ^^^ ^ 
in a suitable neighborhood of Xq. 
Proof of Lemma 13. II 

Let X G Sn- From (3.6), we have for all f = (fi, • • • ,fAr) e "H^ and for all 

ie{i,--- ,N} 

{%iUr) = (f„G(r-.))as+E(fj-'G(x.-x,+r-.))aB- 

Since the first term does not depend on x, we only need to prove the analyticity of 
* : R3\2B -^ C{'H-^l'\dB),'H^/^{dB)), 

f.^g(r):=(f,G(z + r-.))aB 



which follows from the analyticity of G in R^ \ {0} since this formula does not involve 
the singularity of the Green kernel G. D 

We have the following consequence. 

Corollary 1. Since 7^ is an isomorphism, for every x G Sjq, the mapping 
X I — > TxT^ "is also analytic from Sjy to CiJiJ^ i'^n )• 
Proof 
For X G Sm and for |y| small enough, writing the Taylor series of 7i+y as 

Ti+y = 7i + ^7;,5(y), 
where 7i.g(y) denotes the g— homogeneous term in y of 7i+y, we have 

<?>! / J [p>0 \g>l 

where the last expression may be rearranged as a converging Taylor series. D 

3.2. Approximation for large distances. We will also use the asymptotic 
behavior of T~^ as (5(x) :— miuj^j |xj — Xj| tends to oo. For i,j e {1, • • • , N}, we 
introduce the notation 

^i,j • — ^i ^j 1 '^i^j ■ — l"^^ijl an CI i,j ■ — 5 

and, for every f e n^^/'^ (83,11^), 

(rf)(r)=/ G{r-(j)i{a)da, \fr £ dB . 

JdB 

— 1/2 

Proposition 3.2. For every f e Hj^ , 

(Txf). = Tf, + -^^— (Id + e,,,®e,,,)(f„Id)35 + 0(r2)||f||^_,/.. (3.7) 
oiTTj ri j N 

For every g G T-LJ^ , x G 5jv and y G (R'^)^, we have 

(T^-'s). = r-'|g.-^5:^(id + e,,®e,,)(r-^g„id)^^l+o(5-2)llg|l„v. 

(3.8) 

Proof 

Let X G Sn and f G "H^ . From (|3.6[) we have for i = 1, • • • , TV and r G dB, 



(TJUr) - rf.(r) + J2 (fj^ G(x,., + r - .)>aB ■ (3-9) 



Using the expansion 



G(x,,,+z) = G(x,j) + 0(5-2)^ 
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which is vahd uniformly for bounded z, we get 



Subs tituting into (3.9), and using (3.2), we obtain (3.7). For (3.8) we remark that 



(3.7) is the asymptotic expansion in x of the analytic operator 7i at infinity. Since T 
is invertible, we obtain the expansion (3.8) by inverting (3.7). D 



In the sequel, we will use furthermore the fact that the spheres are non-deformable 
and may thus only move following a rigid body motion. In this case, the velocity of 
each point r of the z— th sphere is given by 



Vj(r) = Ui + 0;^ X r. 



(3.10) 



Using (3.8) we find the expansion of the total force that the system applies to the 
fluid 



N 



F = 2_\ 67r77aUi 



n 2 ^ 



J2J2^^JU,+OiS-')\\u\ 



i=l 



«=1 JT^i 



where Sy is the matrix defined by 



(3.11) 



(Id + e 



J J >^ '=^i,j J I 



' « J 



while the total torque (with respect to the origin O) is given by 

N 

T = 67r7/a^Xj x u, + 0(l)||u|| . 



(3.12) 



(3.13) 



4. The swimming problem as an affine control problem without drift. 

In this section we show that self-propulsion of the swimmers makes the system that 
describes their dynamics a linear control system. 

Lemma 4.1. The system that describes the dynamics of low Re swimmers can be 
written in the form 



i(l)^tm.P)i. 



(41) 



where the vectorfields {Fi)i<ci<M o-re analytic in (^,p). It is therefore an affine control 
problem without drift where the controls are the rate of shape changes. 

The rest of this section is devoted to the proof of this lemma. This can be seen 
as a generalization of the result obtained in 3 for the three sphere swimmer of Najafi 
and Golestanian. We show in particular how self-propulsion gives rise to the specific 



form (4.1 ) of the equation, describing the dynamics of all these swimmers of Section 

m 

10 



4.1. Self-propulsion. The question we want to address is whetlicr it is possible 
to control the state of the system (i.e. both ^ and p) using as controls only the rate 
of shape changes ^. For this aim, we need to understand the way p varies when one 
changes ^. This is done assuming that the swimmer is self-propelled, and that the 
swimmer's inertia is negligible, which implies that the total viscous force and torque 
exerted by the surrounding fluid on the swimmer must vanish. As we shall see, using 
these conditions, p is uniquely and linearly determined by ^. The condition that total 
viscous force and torque vanish is written as 

N 

E/ '7;7l,p)V,(^,P,r)da(r) = 0, (4.2) 

4 = 1 "'9B 



and 

N 



^ / X,(^,p, r) X T-l^)Mi,V, r) da{v) = , (4.3) 

,_i J dB 



which lead, using (2.1), to a set of linear equations which link p to ^. In all the 
examples we have in mind (see below), this allows us to solve p uniquely and linearly 
in terms of ^ 

M 

;^ = EW'(e,p)6- (4.4) 

1=1 

In the preceding equation {'Vi{^,p))i<i<M are M vector fields defined on TV, the 
tangent bundle of V . 

The state of the system (^,p) thus follows the ODE 



l)=EF.(e,P)6, (4.5) 



M 
^ i—l 

where now for all 1 < i < Af, F,;(^,p) = f ,,, A x I is a vector field defined on 
~ ~ V Wj(C,p) J 

TS X TV (here e^ is the i-th vector of the canonical basis of R^^). 

4.2. The three sphere swimmer of Najafi and Golestanian (3S). Due 

to axial symmetry, the only non trivial equation in (4.2) is the one concerning the 
component of the total viscous force along the axis of symmetry. This reads 



E 

fe=i 



^Afc(C,p)efc + Ao(^,p)p = 0, 



where 



MO ^ E/ e,.(r-Jfi,t^,t^))* 



3 



- / ei-(7'xT?.p)(«^i'ei,ei))^dr, 
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and, similarly, 



IdB 

while the coefficient of p is simply given by 

3 „ 
MO = y] / ei • (7;7l )(ei,ei,ei)) dr 

> 0. 

Positivity of the viscous drag Ao arising from a rigid translation at unit speed is 
a consequence of the positive-definiteness of the Oseen resistance matrix, see [T7] . 
The coefficients above are independent of p because of translational invariance of the 
problem. From the analyticity of the coefficients and 7^ in x, we have that 

2 
k = l 

where Vk{S,) = — , are analytic functions of ^1,^2- The complete system reads 

Ao(0 

-|(M=«i(i)( \+a,{t)[ 1 I (4.6) 

where aiit) = S,i{t) and a2(i) = ^2(0 ^^^ the control functions. This is under the 
form (|4.1|). 



4.3. The three sphere swimmer moving in a plane (3SP). Due to the 

symmetries of the system, the swimmer stays in the horizontal plane. Therefore, the 
third component of the total force F^ vanishes, while only the third component of the 



total torque T3 might not vanish. The constraints of self-propulsion (4.2 1 and (4.3) - 
more precisely, the vanishing of the components Fi, and F2 of the force, and of the 
third component T3 of the torque - lead to a linear system that can be written as 



Res(C,0)( I )+J2v^it^)^^=0, 



where neither Res, nor V^ depend on c because of translational invariance of the 
Stokes problem. 

The complete system reads 

^ ^ / Id 



^1^1 -Res-i(e,0)V(C,0^ '"(')' (^-^^ 



where V(^, 6) is the 3x3 matrix whose columns are Vi(^, 9), and a{t) — £, G R'^ are 



the controls. This is again under the form (4.1 1 
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4.4. The four sphere swimmer moving in space (4S). Constraints of self 
propulsion (4.2) and (4.3) now become 



K(^,7^) c^(^,7^) 



where the 6x6 matrix 



Ev,:(C,7^)6 = o, 



Res(^,7^) 



c(c,7^) ni^n) 



(4.8) 



known as the resistance matrix jl7) is symmetric positive definite, and depends ana- 
lytically on {^,TZ), just as the vectors Vi{^,TZ) do. 
In this case, the complete system reads 




Id 



4,4 



-Res-^(^,7^)V(^,7^) 



a{t), 



(4.9) 



where V(^,7?.) is the 6x4 matrix whose columns are V.j(^,7?.), and a{t) e R* are 
the controls. Using TZ ~ skew(a;)7?. enables us to write the system under the more 



classical form (4.1 1 



c I =^F,(e,7^)a,(t), 

7^ / i=i 



(4.10) 



where ai{t) — ^i{t) are the control functions. 



5.1. Geometric control theory. When we use (<^i)i<i<A/ as control variables, 
systems like (4.5) belong to the class of linear control systems with analytic coefhcients 
of the form 






M 



Y.a,{t)F,{X) 



(5.1) 



where X is a point of a d— dimensional manifold Ai, and (-Fi)i<i<M are p vector fields 
defined on TA4, the tangent bundle of Ai. For such systems, the tools of geometric 
control theory can be applied. In our case, A4 and the vector fields Fi are furthermore 
analytic, which enables us to use stronger results. The theory for such systems may 
be found in [HITQ], for instance. 

In (5.1 ), the control which governs the system is the vector (ai, • • • , ap), and the 



question of controllability of the system can be stated as follows: 

Question 1. For any pair of states {Xq,Xi) S M^ , are there M functions 
ai : [0, T] —> R such that solving (5.1) starting from X{Q) = Xq leads to X{T) = Xi ? 



A classical result by Rashevsky and Chow states that the system (5.1) is locally 
controllable near Xq e M (which means that question Il| possesses a positive answer 
for any Xi in a suitable neighborhood of Xq) with piecewise constant controls ai that 
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are moreover continuous from the right, if the Lie algebra generated by {Fi, ■ ■ ■ , Fm) 
is of full rank at X^ 

dim Lie{Fi,- ■■ ,FM){Xo)=d. (5.2) 

Here, the Lie algebra Lie{Fi, ■ ■ ■ ,Fi\i){Xo) is the algebra obtained from the vector 
fields {Fi, • • • , Fm) by successive applications of the Lie bracket operation defined as 

[F, G]{X) = {F ■ V)G(X) - (G • V)F(X) . 

The rigorous proof of Rashevsky-Chow's theorem may be found in [TJ[Tn] but we would 
like to give here a hint of why such a result holds. Given a vector field F : A4 — > TM, 
we call exp{tF){Xo) the solution at time t of 

^=^(^)' (5.3) 

X{0) - Xo . 

The main idea is that one can reach from Xq points in the direction g — ^^ f3iFi{Xo) 
by using (5.3 1 with F = ^- PiFi since we have, to first order, for small |i| << 1, 

exp(tF)(Xo) = Xo+ tFiXo) + 0(1^) =Xo + tg + 0{t^) . 

This formula shows that the global displacement of the solution of the dynamical 
system is proportional to time for these directions. 

A more subtle move allows one to reach points in the direction [Fi, Fj]. Namely, 
we now need to consider 

exp(-iF,) o exp{~tFj) o exp(ii^,) o exp{tFj){Xo) = Xa + f[F,,Fj]{Xo) + 0{t^) , 

which also shows that in time i, one can reach a displacement in the desired direction 
which is proportional to t^ . Lie bracket directions of higher order are also attainable 
at the price of higher order expressions in t, leading to smaller displacements. If the 
Lie algebra has a full rank, every direction in M is attainable and the system is locally 
controllable. 

We will focus on systems for which the full rank condition (5.2 1 is satisfied with 
first order Lie brackets only, and a necessary condition for 

dimSpan(Fi,--- , Fm, [Fi,F2], [Fi, F3], • • • ,[Fm^i,Fm]){X^) = d, 

is that 

,^ M(M-l) 

M + — ^— '- > d . 

Our three systems are designed in such a way that they all satisfy MH ~ d, 

and local controllability near X^ follows from verifying that 

det (Fi, • • • , Fm. [-^1,^2], [Fi, F3], • • • , [F^-i, Fm])(Xo) ^ . (5.4) 

The remaining of this section consists in the detailed analysis of our three model 
swimmers. In each case, we show that (5.4) holds at one point, proving local control- 
lability at one point (^0:Po)- In order to obtain global controllability, and theorem 
|2.1[ we need a special argument which is the aim of Section |5.5[ The regularity of the 
shape function ^(t) follows from the fact that the controls are here the rate of shape 
changes. 

14 



5.2. The 3S swimmer. For the equation (4.6), the controhabihty condition 

dVi 



(5.4) becomes 



det(Fi,F2,[i^i,i^2])(C) 



dV2 



(0 



d^ 



(0^0, 



(5.5) 



which is enough to verify for one configuration. For this purpose, and in order to 
further simphfy the computation, we make the observation that, due to symmetry, 

Ai(ei,6) = -A2(6,ei), while Ao(a,6) = Ao(6,Ci)- 



Therefore, if we consider a symmetric shape ^ = (CiC)i showing that condition (5.5) 
holds at f is equivalent to proving that 



96 



(O^O. 



Using the expansion for large distances of the total force (3.11 ) with Ui = (jj — Ci)gi, 
U2 = pei, and U3 = (p + 6)ei, leads to 



Ao(6,6) == STrayy f 3- ^ 



Ai(6,6) = STra?/ f 1 



3a 

y 



1 1 

1 



6+6 
A2(6,6) = 6™,(-l + -(^ + — 



+ o(r'), 



Therefore 



96 



CO-'-^ + oic^) 



which does not vanish for ( sufhciently large. This proves the global controllability of 
the system. 



5.3. The 3SP swimmer. We check the controllability condition (5.4) for the 
equation (4.7) at ^ = and at a symmetric shape ^ = (6 (, C)^ with C e R sufficiently 
large. 

Calling Fj = — Res^^(6 ^)V,j(6 9), the special form of our system (in particular 
the fact that the coefficients do not depend on the position variable c) allows us to 
rewrite the sufficient condition (5.4) for controllability at {£^,c,6) as 



A = det 



9Fi 9r2 aF2 dFs dF3 ^Fi 

a6"96'a6"a6'56~96 



(Ot^O. 



(5.6) 



Notice that the absence of d/d6 terms come from the fact that {Fi)^ = in a sym- 

dFk - 
metric configuration. We now turn to the computation of -rr^iO ioi k y^ I or, at 

least, of its asymptotic expansion in terms of negative powers of 6 
From the definition of F^ , one has 



— — = -Res (6 9) -wy- + ^^^] 
96 V 9£,i d^i 
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(5.7) 



and we need to compute all the factors of this expression at (CiO). 

The first step consists in computing an expansion for Res(^, 0). We get from 



(3.11) and (3.131, and taking N ^ 3, 



T, /c-^^ « /3Id + 0(C-i) 0(1 

Res(e, 0) == 67ra?7 ( ^ 3^ + C 



3C + 0(C) 



and 



1 



1 , -Id + OiC') 
Res-^(e,0)=^^ 3 



3e 



0{C 



(5.8) 



(5.9) 



Differentiating Res with respect to ^; gives (notice that from Section 3.1 Res is 
analytic in ^ and we may therefore differentiate its expansion term by term) 



^(e,0) = 67rar;( "2 
9^1 V (t,xe3)*+0(C-i) 



«4(E.^,S,,)+0(C-3) t,xe3 + 0(C- 

0(C) 



(5.10) 



where S^, originally given by (3.12), is now seen as its restriction to the horizontal 



plane, and the vector t/ x 63 as a horizontal vector. 
Similarly 



gives 



Fk{lO) = -, 



1 / tk + OiC') 







(5.11) 



(5.12) 



(here we have used the fact that, for a symmetric configuration, the displacement of 
one sphere induces no torque) and 






(C,0)= -971776 



^t, + o(c-^) 

o(c-i) 



(5.13) 



From these expressions, and using (5.7), and the identity 

no . 1 

'" = -Y {-{^ii ■ t^Id + t; (g) e^; + e^; (g) t/ - 3(ei/ • ti)eii (g) en) 



one gets (after some calculations 
dF 



fc..-.^ dFi^-^ /#(t/-tfc) 



a&«-»'-i;«>°'n-., 



1 t xt ^ I+OK-'), 

— 7^721-1 X tfc • e3 



(5.14) 



The determinant A given by (|5.6|) is thus equal to 

A 



< 



'^ dct 



t2 — ti t3 — t2 ti — ts \ r)(/—T) 
t2 X ti • eg t3 X t2 • eg ti X tg • eg / 



which does not vanish for C large enough. 

We emphasize again that the global controllability result proven here means that, 
by using suitable controls, the swimmer is capable of moving anywhere in the hori- 
zontal plane and of rotating around the vertical axis by any angle. 
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5.4. The 4 S swi mmer. As before, we check the controUabUity condition (5.4) 
for the equation (4.101 at a symmetric shape ^ = ((, (, (, C,), with ^ sufficiently large. 
Also, we use 7?, = Id in the verification of the Lie bracket condition (5.4 1. 

Notice that for the symmetric configuration ^, the vectors F^ defined by ( 4.9|4.10 1 , 
have no torque components, and that F^ does not depend on c because of the trans- 
lational in vari ance of the problem. Taking into account these facts, the sufficient 
condition (5.4) for controllability at (i^, 0, Id) simplifies to 



dim Span 



5Fi 
56 



dF2 



9F3 



dFi 

96 



(6 Id) =6. 



(5.15) 



Notice that in the preceding equation, each vector is a tangent vector to the group 



of 3 dimensional rigid transformations, SE{3) 



-1IJ3 



S0(3). In ord er to simplify the 



computation, it is convenient to call (p^ = —Iles~^{S^,TZ)Yi (see (4.9)), and work 
in translations and angular velocities. A simple change of coordinates shows that 
condition (5.15) simply rewrites 

'901 



det 



96 



96'' 



903 

96 



904 

96 



(6id)^o. 



(5.16) 



We now turn to the computation of -7^t-{^, Id) for k ^ I, or at least of its asymptotic 

expansion in terms of negative powers of 6 
The definition of 0^, still gives 



90fc 
96 



= -Res"l(67^) 



/9Vfe aRes 



V96 



Oil 



0A 



(5.17) 



and we need to compute all the factors of this expression at (6 Id). 

The first step consists in computing an expansion of Res(6Id). We have from 
(SUl) and dslsl) 



K == 67rar/iVId 



97rr/a^ 



N 



5:^s, + o(r2), 



i=l j^i 



N 



C ^ dnar] ^ 6 skew(tj) + 0(C°) , with skew(t)x = t x x , Vx e R^ 

N 

n = dnar] ^ £_f{Id - t^ (g) t,) + 0(C) • 



i=l 



This gives (here A^ = 4, J2i=i ^i 



OandEti(Id-t.®t,) = |ld) 



Res(6 Id) = dirar] 



and 



Res-i(6Id) 



Gnari 



4Id + 0(C-i) 

iid + o(c-i) 

o(r^) 
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se 



o(C-i) 
Id + 0(0 



(5.18) 



8C^ 



, Id + o(r') 



(5.19) 



Differentiating Res with respect to ^i gives (notice that from Section 3.1 Res is 
analytic in ^ and we may therefore differentiate its expansion term by term) 



^(e-,id) ^ e™, -i«4 i^^^' ^-) + ^(^"') -k«-(*') + ^('^"^^ . 

^^' V skew(to + o(r ') o(c) y 

(5.20) 



Similarly 



leads to 



V. = 



eTTTyatfc - 97r?7a2 X;,^fc Sjfetfc + 0(C ^) \ , . 

o(i) j ^^■^^> 



0.(C-Id)^4(^^- + ^(^-^)) (5.22) 

(here we have used the fact that for a symmetric configuration, the displacement of 
one sphere induces no torque) and 

"6 \ o{c') 



From these expressions, and using (5.17), one gets after a tedious but straightfor- 
ward computation 

^'^fc (^- Id) - |^(e,Id) = ( 6^^*' " ^'^O + OiC') ■ (5.24) 



56 ' d^k ' \ ji^ti X t 



lec^ 



Cfc 



The determinant A given by (5.16) is thus equal to 

A = C''(^^\dct(l^-l' t3-ti ... t4-t3 \ 13) 

1^1024^/2^ V^axti t3 X ti ... t^ x t^ J 

which does not vanish for C, large enough. 

5.5. Prom local to global controllability. In all the three model swimmers 
that we proposed, we have found a point {S,o,po) (z S x V at which the sufficient 
condition for local controllability 

det(Fi,... ,Fa,,[Fi,F2],... ,[Fm-i,Fm]){^o,Po)^0, (5.25) 

holds, which proves, denoting by T the family of analytic vector-fields (Fi, . . . , Fm), 
that 

dim Lie(^„,p„)(^) = M+ ^"^^^^ ~ ^^ = d. 

In order to show the global controllability of the problem, we show the following 
lemma 

Lemma 5.1. For every {^,p) <E S xV, 

dim Lie(^p){T) = d. (5.26) 
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Proof 

Let { £,,p) € S X v. From Hermann-Nagano Theorem (see [12], P- 48), we know 
that (5.26) holds at every point {£,,p) in the orbit Ojf(^OjPo) of J-" passing through 
{£,o,Pq)- Now, due to the specific form of the (Fi)i<i<M, there exists q £ V such that 
(C,g) e Ojr{Co,Po)- Indeed, 

/A/ \ 



i^,q)=cxp ^(^-^o),:F, i^o,Po) ■ 



Eventually, because of the invariance of the problem with respect to the orientation 
of the swimmer, one has 

dim Lie^^p-j{J-') = dim Liet^^q){F) = dim Liei^^g^p^){F) — d. 

This ends the proof of Lemma |5.1| and thus of Theorem |2.1[ D 



6. Optimal swimming. A possible notion of optimal swimming consists in 
minimizing the energy dissipated while trying to reach a given position px at time T 
from an initial position po at time 0. In this sense, the optimal stroke is the one that 
minimizes 

l-T N N 
£i^^P)--^ EE/ MLp,r)-T~^lp^v,{^,p,r)da{r)dt (6.1) 

"'0 i^i j^i JdB 



where v satisfies the self-propulsion conditions expressed by Equations (4.2 ) and (4.3 1, 
subject to the constraint that p(T) is equal to px- 

The special structure of our problem allows us to express the dissipated energy 
solely in terms of the shape variables £,, by virtue of equation (4.4). We start by 
expressing the velocity of the tti— th ball due to a unit shape change in the i-th 
component of ^ as 

wr(C,p,r) := V^X„.(e,p,r) • e, + VpX„,(^,p,r) • (V(e,p) • e,) . (6.2) 

This allows us to simplify the energy of the system to an expression depending 
on the rate of shape and position changes ^: 



„T M M 

£i^^P)--= EE^^^^^(^'P)0'^i> (6-3) 



where the metric Q is defined by 

N 

G^jiLp) ■■= E / wr(^,P,r).r^7i (w](e,p,r),... ,wf(^,p,r))™da(r). (6.4) 

Since dissipation is invariant with respect to the current position p, it is possible 
to show that Qij{i,p) = QijH, 0) := Gij{C) for any position p. 

Optimal swimming is then given by the solution {£,^p) : [0, T] — > S x V oi 

min / i-g{0-idt^:£{0 (6.5a) 

subject to C(^,p,i)=0 Vie[0,r] (6.5b) 

a<m<^u Vie[o,r], (6.5c) 
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where the contramts C{S,,p,t) — are given exphcitly by 

V<e[0,T], 



Vim, Pit)) -m- 


-Pit) =0 


m^^T) 


= 


P{0)-Po 


= 


p{T) - PT 


= 


m - eo 


= 



(6.6) 



Sonic of the constraints in (6.6) can be relaxed. For example, it is possible to fix 



only the periodicity of the shape, but not its actual initial and final values (so that 



the last constraint in (6.6) disappears), or to fix only some of the components of p(T), 



leaving the others free to be optimized. 



7. Numerical Approximation. Our strategy for the solution to problem (6.5) 
is to discretize in time both ^(i) and p{t) using cubic splines, and rewrite everything as 
a finite dimensional constrained optimization problem, defined only on the coefficients 
of the spline approximation of ^ and p. 

A C++ code for boundary integral methods, which has been developed using the 
deal. II library (see 8J and [T), solves the Dirichlet-to-Neumann map T^^, and is 
then used to construct Q{£,) and V(^,p). 

The finite element method in the spline space is then used to produce a finite 



dimensional approximation to (6.5 ), which is solved using the reduced space successive 
quadratic programming algorithm (rSQP) provided with the Moocho package of the 
Tri linos C++ library (TS]. 

In this section we touch only briefly on the numerical implementation of the 
discrete Dirichlet-to-Neumann map for a collection of spheres, and we refer to [2] for 
more details on the rSQP technique used for the search of constrained minima. 

In principle, standard surface meshing techniques could be employed to integrate 
on the surface of our swimmer, but they lack the speed and accuracy that specialized 
methods can offer when integrating on such special domains like collections of spheres. 

The method we have chosen for the integration on the spheres is presented in [l6] 
for use with smooth functions and experimental data. It uses an oblique array of 
integration sampling points based on a chosen pair of successive Fibonacci numbers. 

The pattern which is obtained on each sphere has a familiar appearance of inter- 
secting spirals, avoiding the local anisotropy of a conventional latitude-longitude array. 
We modifled the original method to take explicitly into account multiple spheres and 
the singularity of the Stokeslet (see [28] for a more detailed introduction on boundary 
integral methods for viscous flow). 

The discrete versions u/j and f/; of the continuous velocity and force densities 
u and f are obtained by sampling the original functions at the quadrature points, 
translated and replicated on each sphere according to ^ and p, to obtain Ny scalars 
that represent Nv/3 vectors applied to the quadrature points sampled on dfl. 

Integration on dQ of the function u^ (x) is then effectively approximated by 



/ 

J on 



Nv/3. 

u(x) dcr(x) - Y^ u(q*)o 



where the weights oj' and quadrature points q' are derived according to 
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The construction of a discrete Dirichlet-to-Neumann map for our swimmers fol- 
lows the usual collocation method (see, for example, |23]), 

u(q^)= / G{y-q')i{y)da{y) 
Jdil 
Nv/3 _ (7.1) 

- Y, G(q-'' - q*)f (q-'' V + G"f (q*y i = 1, . . . , 7Vy/3. 



A specialized calculation of the singular integral is performed in Equation (7.1) 
when i is equal to j, by deriving G according to the (known) exact solution of the 
flow around a sphere. 



Equation (7.1) can be expressed more compactly as: 

{h=A-'uh, (7.2) 

where the matrix A is obtained by writing explicitly the various components of Equa- 



tion (7.1), and we have used the notation Uh and fh to indicate the Ny dimensional 
vector that samples the functions u(x) and f (x) at the quadrature points q. 

The construction of the discrete Dirichlet-to-Neumann map allows us to compute 
explicitly, given a pair of shape and position variables {^,p), the numerical approx- 
imations of GhiO ^^'^ ^ft,(CjP)- These can in turn be used to solve the constrained 
minimization problem as in |2], i.e., by using the implementation of the rSQP method 
included in the Moocho package of the Tri linos library [5]. 

8. Numerical Results. In this section we analyze some of the optimal strokes 
that can be obtained with our software for the case of the three-sphere swimmer in 



the plane [23] (3SP), depicted in Figure 2.2 While this is one of simplest possible 



extension of the three-sphere swimmer by Najafi and Golestanian [27], it adds already 
a high degree of complexity to the numerical approximation scheme. 

The added complexity comes both from the presence of three additional time 
dependent variables (one shape and two positions) , and from the fact that the vector 
field that relates positional changes to shape changes {p — V(^,p) • f) depends also 
(though in an explicit way) on the orientation 9 = p^ oi the swimmer. 

In the software used to compute the optimal strokes for axisymmctric swim- 
mers [2], there is only one position variable pi{t) — c{t), which can be obtained by 
post-processing the vector field V(^) (which, in the axisymmctric case, is independent 
on the position p) and ^ itself. 

In the present case of a planar swimmer, this is no longer possible, and the 
optimization software has to go from optimizing over the two shape variables ^i(t) 
and £,2{t) to optimizing over the entire shape-position space, composed of the six time 
dependent variables ^(i) = i^i{t),£,2{t),^s{t)) and p{t) = {cx{t),Cy{t),6{t)). 

In the experiments that follow, we use water at room temperature (20° C) as the 
surrounding fluid, and we express lengths in millimeters (mm), time in seconds (s) and 
weight in milligrams (mg). Using these units, the viscosity of water is approximately 
one (ImPas = Img mm~^ s~^), and the energy is expressed in pico- Joules (IpJ = 
Im.gm.m'^ s^^ = 10^^^ J). 

The radius of each ball of the swimmer is .05mm, and additional constraints 
are added to the optimization software to ensure that each component of the shape 
variable ^ stays in the interval [.Imm, .7mm\, so that the balls never touch each other 
and that they don't separate too much. 
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We present some of the computed optimal strokes that iUustrate the richness 
of the problem. In all the examples that follow, we ask the swimmer to go from 
p(0) = (0mm, 0mm, 9o) to p{T) = {.Olmm, 0mm,, 0t), for various combinations of ^o, 
0T andC(O) =^{T). 

Since the swimmer is 27r/3 periodic, it is sufficient to study the optimal strokes 
varying 9q between and 2tt/3 while keeping c^ and Cy constant to explore all possible 
optimal swimming strokes. All other combinations can be obtained by rotations and 
permutations between the components of the optimal solution (^,p). 

Moreover, since rotating the swimmer by n is equivalent to asking it to swim 
in the opposite direction, reversibility of Stokes flow allows us to further reduce the 
study of the optimal strokes to the interval 9q S [0,7r/3]. 



In particular, if the pair ^(t) and p{t) = (c(t), 9{t)) is an optimal solution of (6.5), 
then we have that all the following permutations, symmetries and rotations of {^,p) 
are also optimal swimming strokes, but with different initial and final positions: 



m 


c{t) + k,9{t) 


VfceR^ 


(8.1a) 


m 


7^,(7)c(f),e(t)+7 


V7 e [0, 2tt] 


(8.1b) 


Ut),Ut),Ut) 


n,{2Tr/3)c{t),9{t) 




(8.1c) 
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Given a family of optimal solutions with different starting angles ^o G [0,7r/3], 



then properties (8.1 1 can be used to derive the optimal solutions in all of the interval 
[0,27r]. 

The optimal strokes we present have been formatted to show the trajectory of the 
center of the swimmer amplified by a factor 60 with respect to the dimensions of the 
swimmer itself, and the orientation 9 of the swimmer amplified by a factor 10. The 
full configuration of the swimmer is shown at ten equally spaced times between t = 
and i = T = 1. 



Our first optimal stroke (Figure 8.1 ) is obtained by fixing the initial and final rota- 
tion 9q — 9t to be zero and the initial shape of the swimmer to be ^ = (.4, .4, A)mm. 
The total energy consumption to swim by .0\mm in the x direction with a single 
stroke is 0.5 lip J. 

The first important thing to notice is that, although the initial configuration is 
symmetric with respect to the translation axis, the stroke is not. This implies that 
the optimal solution is not unique: exchanging the off axis shape variables ^2(i) and 
S,'i{t) yields to the same final displacement with the same total energy consumption. 

A configuration which is symmetric with respect the y— axis can be obtained by 
setting ^0 = '"'/S, and also in this case the solution is non-symmetric and non- unique. 



Given the permutation properties (8.1), we deduce that two distinct global so- 



lutions with the same energy dissipation exist for 9q — Ntt/6, for any integer N. 
From the existence of these different global minima we can deduce that at least two 
branches of local minima exist outside those particular points, starting off as contin- 
uous perturbations of the two global solutions. 

Numerical evidence show that at least three stable branches exist, as shown in 



Figure 8.2 where we plot the energy dissipation of all the (local) optimal strokes with 
respect to the starting initial angle ^3(0) = ^o- It is possible to steer the software 
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Fig. 8.1. 3SP swimmer: optimal stroke for Bq = = Ox radians. The initial shape is fixed at 
^(0) = ^(T) = (.4, .4, A)mTn. Energy consumption to swim by .01mm in the x— direction: 0.511pJ. 




^/'■^ 0(0) 



Fig. 8.2. Branches of optimal swimming for 3SP swimmer: dissipated energy VS starting angle. 



towards one or the other local branch by feeding it with initial guesses which are close 
to the desired branch. 



From Figure 8.2 it is also clear that the swimmer prefers to swim along directions 
which are close to parallel to one of the arms: Ntt/S ± 7r/30 are the approximate 
preferred swimming directions. The directions perpendicular to one of the arms {9q = 
7r/6 + Nir/S) are the most expensive ones (although they differ in energy only by 
approximately 2%). 
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«0)=e(r) = (.4..4,.4)r, 



r 




Fig. 8.3. 3SP swimmer: optimal stroke for Bq = 0, 9^ optimized to .058 radians. The initial 
shape is fixed at 5(0) = ^{T) = (.4, .4, A)mm. Energy consumption to swim by .01mm in the 
X— direction: .209pJ. 



From Figure |8.1| it is evident that a lot of energy is wasted in forcing the final 
rotation to be equal to the initial one. If we relax this requirement and allow the 
software to optimize the stroke without forcing the rotation to be periodic in time, then 
the swimmer does not follow an "eight-shaped" path anymore, and the actual target 
displacement is reached with a smaller stroke (Figure 8.3 1, and by dissipating less 
than half the energy (.209pJ instead of .511pJ). The final rotation is approximately 
.058 radians (about 3.3 degrees). 

The main defect associated with this strategy is that all successive strokes will 
be rotated with respect to the initial swimming direction. If the stroke were to be 
repeated, we would observe a rotational drift and an average circular motion along 
a trajectory of radius r — .35mm. To ensure that several successive strokes produce 
an average trajectory along a straight line, one cannot use copies of the same stroke. 
Instead, different optimization problems need to be solved for each stroke, in each of 
which the initial orientation is the one reached at the end of the previous stroke. 

A complementary strategy to reduce the energy of a single stroke without gener- 
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t ?(0) = ?(T) = (.15, .15. .19)mm 




Fig. 8.4. 3SP swimmer: optimal stroke for d^ = = dj- . The initial shape is optimized to 
^(0) = i{T) = (.15, .15, .19)mm, Energy consumption to swim by .01mm, in the x— direction: .221pJ. 



ating a rotated final configuration consists in optimizing also on the starting shape, 
rather than on the final angle. Figure |8.4| shows this strategy, where the final angle 
is fixed to be equal to the initial angle. From the figure it is clear that the more com- 
plex hydrodynamic interaction due to the closeness of the spheres can be effectively 
exploited by the optimization software to reduce the dissipated energy. 

The optimal initial shape is found to be (.15, .15, .19)mm, a non-symmetric shape. 
The total energy dissipation with this stroke strategy is about .221pJ. While this is 



already less than half of the energy consumed by our first optimal stroke (Figure 8.1 1 
it is still about 10% more expensive than the second attempt (Figure 



.3) 



The combination of the two strategies, in which both the initial shape and the 
final angle are left free for the software to optimize, is shown in Figure [83) In this case 
the total energy dissipation drops down to .lOOpJ. The initial shape is optimized to 
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f(0) ={(r) = (.23, .35, .2)mm 




Fig. 8.5. 3SP swimmer: optim,al stroke for Bq = 0, 8t optimized to .057 radians. The initial 
shape is optimized to ^(0) = £,{T) = (.23, .35, .2)mm. Energy consumption to swim by .01mm in the 
X— direction: .lOOpJ. 



the value (.23, .35, .2)mm and the final angle is 0.057 radians, or about 3.28 degrees. 
In this case the optimal stroke hits the boundary of the region of allowed shapes. 
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